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I. INTRODUCTION 


A. BACKGROUND 

Extremely low frequency (ELF) geoelectric noise has a 
highly impulsive nature which makes it very difficult to 
remove by filtering. The main contributors to the 
impulsiveness of the ELF spectrum are lightning discharges. 
The lowest frequencies of the ELF range (1 to 60Hz) have such 
a low attenuation rate as they propagate, that lighting 
discharges all over the world must be considered when studying 
the noise spectrum at a single location. ELF signals have 
been studied for such applications as power transfer, 
communications, geophysical surveys, and so forth. The amount 
of information on the subject is immense and increasing. How 
to account for the impulsive background electromagnetic fields 
inevitably present in any system operating in this frequency 
range remains a serious problem that any system designer must 
answer. 

The background noise of the lower ELF frequency range can 
be thought of as consisting of both local and global 
phenomena. The global component is comprised of signals that 
are relatively coherent over extended distances. A detailed 
recording of data in this frequency range 


...1S characterized by bursts of a few cycles duration 
which may be ten or more times the mean amplitude between 


bursts. This bursty character 1s apparently coherent over 
MOSt, Tf NOt all, ot Ee Earrm (res 


The local phenomenon is comprised of those components that are 
not coherent over extended distances and are assumed to be 
products of the local environment of the detector. A set of 
parallel detectors separated by some distance can be 
coherently combined to remove the global portion of the 


background noise and thus lower the noise floor. 


B. OBJECTIVE 

The Applied Physics Lab (APL) at Johns Hopkins University 
has for years studied and made measurements of the background 
geoelectric noise on widely separated pairs of electrodes. 
The purpose of this thesis is to describe techniques of 
adaptive filtering that the author applied to data collected 
by APL, in order to determine the amount of noise 
cancellation. The techniques are implemented using software 
developed during an extended visit to the Applied Physics Lab 
bY sbhe Tauiehiow: In a related effort, software will also be 
presented implementing appropriate detection methodology for 


the detection of narrowband signals embedded in the noise. 


C. ORGANIZATION 

The thesis begins with a brief discussion of the 
constituents of this ELF noise to demonstrate the difficulty 
in modeling of the spectrum and to provide supporting 


background information. The basic concepts of the adaptive 


filtering will be presented as a lead in to the adaptive 
algorithm to be used to cancel the noise. These introductory 
sections will be followed by a discussion of stability and 
convergence concerns for the adaptive algorithm in this type 
of data environment. The detection discussion is then 
presented as a continuation of the processing required for the 
detection of narrowband signals. The adaptive algorithm is 
then used in the detector scheme to demonstrate performance of 
the developed software. A final section will summarize the 
work performed, draw relevant conclusions about the work and 
propose possible future studies that could be derived from 


this work. 


Ii. THE ELF BACKGROUND 


Many studies have been conducted in an attempt to develop 
an understanding of what actually constitutes the natural and 
man-made background noise spectrum in the one to sixty hertz 
range. The sources of the energy in this frequency range can 
be primarily attributed to atmospheric disturbances and man- 
made power transfer systems. A major problem with these 
sources of noise is that they are highly irregular in both 
magnitude and duration of the noise signal. Attempts to 
classify and assign specific characteristics to these sources 
have proven difficult. A brief history of the research 
conducted in these areas is provided in the following two 


subsections. 


A. ATMOSPHERIC NOISE 

The source of the majority of noise in the one to sixty 
hertz frequency range is lightning discharges. The lower ELF 
range discussed here also contains what is called the Schumann 
Range. The Schumann range takes its name from the early 
theoretical work of the German scientist W. O. Schumann. 
Schumann theorized that the cavity created by the surface of 
the Earth and the Ionosphere has naturally occurring resonant 
frequencies. He took this research further to derive a set of 


harmonic resonant equations for calculating these frequencies. 


The work of Schumann has been carried further with the 
actual measurement of the Schumann resonance frequencies in 
field experiments. ie i too And loou ee H. Lt von. Konig 
substantiated the presence of Schumann resonances by observing 
the waveforms in the output of a narrow band amplifier. The 
experimental presence of Schumann resonances created a large 
quantity of literature on the subject as scientists attempted 
to understand the phenomenon. Once a foundation of the 
propagation properties of the Earth-Ionosphere cavity was 
established, the search for the source of the excitation 
began. [Ref. 2] 

The most prominent excitation source of the Earth- 
Tonosphere cavity is the cloud to ground lightning stroke. 
The number of similar strokes present in any one flash 
observed during a thunderstorm is a random variable that 
normally ranges from two to twelve. Since the number of 
strokes 1s generally too high to measure individual waveforms, 
a more applicable measure is the power spectral 
density[Ref. 3]. The frequency spectrum of this noise was 
resolved in sufficient detail to estimate the quality factor, 
QO, of the earth-ionosphere cavity by M. Balser and C. A. 
Wagner in 1960. The Q of a cavity resonator is a measure of 
the bandwidth of the resonator[Ref. 4]. The source of the 
noise, as theorized by H. Raemer in 1961, is the response of 
the earth-1lonosphere cavity to electrical discharges created 


in thunderstorms all over the world. Raemer attempted to 


model the lonosphere in an attempt to reproduce mathematically 
the measured response of Balser and Wagner. Although Raemer's 
attempt fell short of its goal, it began the process of 
refinement of a working model of the ionosphere that continues 
today. [Ref. 2] 

The early work performed by these scientists brought out 
the difficulty in attempting to model this naturally occurring 
phenomenon due to its randomness. Parameters must be 
estimated by some method in order to model this random 
process. The highly variable nature of these parameters makes 
the task of modeling the noise nearly impossible. Examples of 
these varying parameters are diurnal variations in the 
1onosphere, twenty-four hour variations in the ionosphere, 
seasonal variations in the locations of thunderstorm regions, 
and the randomness of individual lightning strikes. These are 
only a few examples of the parameters that must be considered 


in order to model the atmospheric noise. 


B. MAN-MADE NOISE 

The presence of atmospheric noise 1s not the only concern 
that must be addressed when discussing background noise in the 
one to sixty hertz range. Many components in today's world 
emit unshielded interference that falls in this range. These 


components are the result of either faulty equipment or poor 


design. Examples of this noise source are power lines that 
are used to distribute alternating current through populated 
areas. 

Large ferric objects can create fields around them that 
can be detected. The strength of these fields can be well 
below the previously mentioned sources of noise but can be 
discernable in certain Situations. Motion of these large 
Magnetic objects in the Earth's natural magnetic field is the 
source of a small field. 

Motion of the detector can create an increase in the 
background noise. Any slight motion of the detector can 
create misleading information and in most cases, an increase 
Mieseme noise signal. The background ELF spectrum can be 
thought of as the summation of three orthoganal vectors. 
Relative motion between these three vectors and the detector 
creates a false signal that can raise the noise floor. This 
type of noise can be compensated for by the use of motion 
sensor. [Ref. 5] 

Other sources of man-made noise can be found in the 
detector equipment itself. Improper shielding of cables and 
electronic devices used in the design of the detector and its 
processing components can lead to elevated noise floors. 
Unlike naturally occurring atmospheric noise, these sources 
can be traced and eliminated. The removal of this noise can 


Semiewin two forms, elimination of the source or isolation of 


the source from the detector. The fact that man-made noise 
sources exist must not be overlooked when attempting to 


process signals at the output of the detector. 


IfiiL. ADAPTIVE FILTERS 


A. GENERAL TERMINOLOGY 

The ELF background noise has been described as a highly 
impulsive, non-stationary random process. Some of the causes 
of these variations have been briefly described above. The 
use of traditional filter designs may not be the optimum 
method for reduction of the noise floor for reception of low 
strength signals in the one to fifty hertz range. Since their 
introduction, the use of adaptive filters has shown promise in 
combating exactly this type of scenario. To understand the 
application of adaptive filters in this case, a brief 
discussion of relevant terms must be conducted. 

An adaptive filter is by definition a filter that has the 
Capability to adjust, by itself, a set of design parameters 
that are based on estimated statistical characteristics of 
the signal to be filtered. This process of adjusting the 
filter to the present situation alleviates the need to have a 
priori knowledge of the relevant Signal characteristics. The 
previous chapter described how difficult it 1s to derive a set 
of statistical characteristics for the ELF no1lse environment. 
If @ priori knowledge of the relevant signal characteristics 
were known, an optimum filter could be designed. The goal of 


adaptive filtering is to estimate these statistical parameters 


and to refine the estimation through the use of an algorithm. 
If an appropriate algorithm can be found, then after a 
sufficient number of iterations, the adaptive filter should 
converge to the optimum filter. The use of adaptive filters 
requires the designer to look for the optimum algorithm for 
estimating the relevant signal characteristics instead of 
attempting Eo discover the aC Gite signal 
characterist ress Ret ancl 

Adaptive filters are divided into two basic categories, 
open-loop and closed-loop. An open-loop configuration is 
usually a two stage process. The first stage is used to 
"learn" the statistics of the relevant signal. The results of 
the first stage are then used in the second stage to compute 
the filter parameters required using a nonrecursive algorithm. 
In contrast, the closed-loop configuration uses only one stage 
to develop the filter parameters: the relevant signal 
characteristics are not explicitly estimated but are inferred 
from a recursive algorithm that updates the filter parameters 
directly as each new data point is obtained. The adaptive 
filter gains additional knowledge Gis the Signal 
characteristics during each iteration. The gain in knowledge 
results in an improvement of the filter parameters which are 
adjusted and their performance is used as an input to the next 
iteration of the algorithm. The adaptive filter discussed in 
this thesis falls in the closed loop category. The advantage 


of closed-loop over open-loop is that the requirement for only 
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one stage generally corresponds to a less expensive 
configuration based on its simpler implementation. [Ref 6] 

The description of the filters that will be analyzed in 
this study carries terminology that is pertinent to all 
filters, not just adaptive filters. Filters are implemented 
in either continuous-time or discrete-time which defines the 
form of the relationship between the input and output of the 
filter. The filters that will be studied are of the discrete- 
time version which means that the filter may be described by 
a difference equation. The =Structuremorstne fLlleers to be 
Studied are referred to as finite-impulse response (FIR), 
tapped-delay-line, or transversal filters. This description 
means that the filter's output relies only on the past and 
present values of the input. The advantage of a FIR filter is 
its inherent stability. [Ref 6] 

The algorithm that 1s used to update the filter parameters 
1s generally named after some of its operating 
characteristics. The algorithm used here is called the 
Sequential Regression Algorithm or SER[Ref. 7]. pace 
algorithm uses an iterative approach to approximating the 
weighting factors required for a linear combiner. Figure l 
shows a basic diagram of a linear combiner. 

This system uses simultaneous samples of the input signal 
at time index k, to produce an output signal, y,, by a linear 
operation with weighting factors, w,. A second way to 


physically interpret the linear combiner is to look at it in 


ak 





Figure 1 Basic Linear Combiner 


the form of a transversal filter. Figure 2 shows how this 
system 1s designed: X, corresponds to the last L samples of 


xX at time index k. The weighting factors are dual indexed on 





Figure 2 Basic Transversal Filter 


the time delay, 1, followed by the time index, k. The weight 
factor, w,,, corresponds to the weight factor for the sample 


taken two samples ago from the present sample time, k. 


TZ 


B. THEORY AND DESIGN 
1. Configuration 
The filter to be developed here is a combination of 
the previous two figures. The transversal filter is applied 
EOra set of multiple sensor outputs to form an input vector 
that is N x L data points long, where L is the number of 
separate inputs or sources, and N is the length of the 
transversal filter. There are two reasons for using a filter 
design of this nature. The first is the ability to use more 
than one source of input, allowing for a better spectral 
sampling of the background ELF energy. The second is the 
ability to use the previous N samples thereby smoothing the 
output since it is not dependent on instantaneous values which 
may vary greatly from sample to sample. Both of these ideas 
are beneficial to the overall performance of the filtering 
process. 
The output of the filter, y,, can be expressed in 
terms of the input matrix and the weighting factors matrix. 
To do this, some defining equation must be expressed for the 


Maeteninatrix, X,, 


me. = eek. eee ey (1) 


and weight factors matrix, W,, 


W,=[Wi, Wy, oes Weyl? (2) 


where each sensor's input vector is 


the 


Ay ms exec) XxX (Kay) xX, (K-2) soe Xx, ie (N= Rl (3) 


and weighting factor vector is 


Wri = (Wig Wrceeay Waceeay 6+) Wice-cn-ay)l - (4) 

The first subscript of Equations 3 and 4, labeled [zy 
., L, indicates the source of the sample. The second 
subscript, k, indicates the time of the sample. A total of N 
sample values from each sensor are used for each calculation 
of y,. A physical model for the filter to be studied is shown 
ify “Pigure 37 The scalar output value, y, in Figure eee 


defined as: 


Vx X,W,.= WiX,. (5) 

These basic equations will become the starting point 
for the development of an operating filter system for the 
removal of the ELF background noise. It ‘s »reterable that 
the weighting factors, W,, cause the output, y,, to be an 
estimate of some desired signal, d,. The concept of the 
filter is to use the weighting factors to model the 
differences between signals originating from the same source 
but sensed at various locations. The inputs, xX, in Figure 3, 
represent samples of a signal taken at different locations. 
The desired signal, ad, 1s also a sample of the same 
background source signal. The difference between the desired 
Signal and the output of the filter is expressed as the error 


in the estimate, €,. 
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Ex = ay ~ Vx (6) 
In the system designed here, the estimate error is the 
eueout sought for analysis. If the weight factors are set 
properly, they will theoretically remove all signals that are 


coherent between the different sample locations. The goal of 





Wi ec-(-1 
woe Wi ge-(n-1)) 





Vx 


Figure 3 Multi-Channel Transversal Filter 
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such a system is to remove the correlated background noise 
from a Single sensor thus rhaneeine Signals that are 
particular to that sensor alone. If the estimated value of y, 
can be made to represent the ELF background noise, then the 
estimate error, €,, can be assumed to be due only to the 
Signal present in the desired signal. The resulting output 
from the noise canceler has a much lower noise floor due to 
the removal of the correlated noise. 
2. The Mean Square Error 

The filtering process must have a means by which an 
optimal solution can be defined. The mean-square error will 
be used to develop a defining equation, called the performance 
surface, for the optimal solution. Combining Equations 3 and 
4, a definition of estimated error in terms of the input 
vector, the desired signal, and the weig:.ting factors is 
obtained. The expected value of the squared estimated error, 
E(é€,)*, assuming that €,, d,, and xX, are statistically 


Stationary 1S given as: 


E(e4] = El dZ]+W'E(X,X,JW-2E[d,X;]W (7) 
The variables x, and d, are assumed to be dependent, and in 
this case both contain the same ELF background. The equation 
(ONE the mean square error, MSE, associated with the estimate 
is thus defined in terms of the auto-correlation function of 
the input vector and the cross-correlation function between 


the input vector and the desired signal. 
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The auto-correlation and the cross-correlation 
functions will be represented by R and P respectively. The 
elements of these two matrices are constant values if the 
input vector x, and d, are stationary random variables. A 
mathematical expression for the auto-correlation matrix, R, 


oe 


2 
Xik X1 Xo a Niki 
2 
X11 k X2k oie Xo Ark 
R= £ (X,X;)=E] ° (8) 
2 
X11 k Xrikik oy Xrk 


eimfemetehe cross-correlation matrix, P, is: 


P= Eld,X,] = ElAX, AX, ++ AX? (9) 
Using these two definitions in Equation 7 results in the 
following definition of the mean square error in the k*™ 


estimate of y: 


MSE +J,,=Ele;] =E[dg] + W7RW - 2P7W (10) 
The optimum solution for the selection of the weighting 
factors is the condition where the mean square error is at a 
minimum. 
3. Minimizing The Error 
The minimum mean square error condition can be solved 
for by taking the gradient of the MSE with respect to the 


weighting factors and setting it equal to zero. The use of 


Ly 


the minimum value of the MSE as the optimal value 1s based on 
the assumption that the weighting factors are only designed to 
model that portion of the input vector that 1s correlated with 
the desired signal. If the error is at the minimum value, 
then the filter is removing the maximum amount of correlated 
data between the input and the desired signals. The system 
being evaluated is based on the assumption that only the 
desired signal contains both noise and a signal of interest 
while the input vector 1s composed of a pure noise signal 
alone. Under this assumption, the only correlated data 
between the two input data samples is the noise. The desire 
1S to remove as much noise as possible and thus to remove the 
maximum amount of correlated data between the input vector and 
the desired signal. It is therefore assumed that the filter 
1S working at its optimum value when the weighting factors 
meet the criteria of minimizing the mean square error of the 
GueO UE 
4. The Optimum Weighting Vector 

The optimum weight factor will be represented by the 

symbol W°’. The two conditions that must be met in order for 


the value of J to be minimized are: 


Vu Jes lw-wo = 0 (11) 


and 
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Ow ,Ow, 
These two equations state that when the weight vector is at an 
optimum the performance function is at a local minimum since 


the slope is zero and the curvature is upward. The gradient 


of the performance function is evaluated to be: 


V = V,J=V,Eld/] -2V, (WTP) +V,, (W7RW) 3) 
=-2P+2RW 
Setting the gradient to zero and solving for the optimum 


weighting factors results in: 


Wwe = Rip (14) 


This equation is sometimes referred to as the Weiner-Hopf 


equation where the weight factor is referred to as the Weiner 


weight vector. The input autocorrelation matrix, R, can be 
mieerm@ed 1f it 18 "positive definite," i.e., if: 
V'RV > 0 (15) 


where Vis any non-zero vector. When the matrix Ris “positive 
definite" then the optimum weight factors can be solved for 
directly if the complete statistical properties of the auto- 
correlation and cross-correlation functions are known. The 
statistical properties are not known in this case so further 
development is required to obtain an estimate of the optimum 


weight factors. 


Eg 


The second condition that ensures that the error 
estimate 1S at a minimum is that the partial derivative of the 
gradient with respect to the weighting factors 1s greater than 
Zero. To prove this the we take the partial derivative of 
Equation 3 with respect to W 


Cos 
Ow , Ow, 


= 2R (16) 
This result again shows that the input autocorrelation matrix 
must be "positive definite" in order to solve for the optimum 
Weighting factors. 

The performance surface can now be expressed in terms 
of the minimum error possible, the weighting functions,and the 
autocorrelation of the input vector. The minimum mean square 
error possible can be found for the case when the weighting 
vector 1s at its optimum value. The optimum weighting factor 


1s used to solve for the minimum error by substituting W into 


the minimum error equation. 


E[dz] + [W°] TRW°-2P7wW° 
E[d¢] + [R?2P]7RR'P-2P7R Pp 
E[dz] - P™R'P 

E[d;] - P™w° 


Grin 


I 


(172 


Nf 


The minimum error can now be used as the base line by which 
all other errors are found. The error at any time is the 
minimum error with the addition of uncertainty dined 


components used to derive the minimum error. Thus: 
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J = Ini, + (W-W°) TR(W-W°) (18) 
The expression can be shown to equal that of Equation 10 by 
the following proof. Using the fact that (A-B)'’ in general 


can be expressed as A’-B': 


J = J..,+ [(W°] "RW°+W7RW-W’RW? - (W°] "RW deel) 
Since all terms are scalars, the transpose is equal to the 


scalar value, thus the last two terms are equal. Substituting 


for the minimum error and the optimum weighting vector: 


J 


E[{dZ] -P7R P+ [W°] TRW°+W7RW-2W’RW° 
E(d/] -P’R'?P+P7R' RR P+W’RW-2W’RR?P 
= E[dZ] +W’RW-2W’P 

E[d/Z] +W7’RW-2P'W 


(20) 


This result validates Equation 18 by showing that it 
equivalent to Equation 10. An important note about this 
expression of the performance surface is that it 1s quadratic 
with respect to the weighting vector. 

5. Iterative Calculation of W, 

The optimum weight factors can be arrived at 
iteratively by use of a gradient search of the performance 
function. What this means is that after each estimate of the 
optimum weight vector, W,, the slope or gradient of the 
Bem~rormance function, in the direction of the optimum 
weighting vector, is used to derive the next weighting vector, 


W,,,. This method ensures that each estimation of the optimum 


2a 


weight factors is at least as close or closer to the optimum 
value than the previous estimate. Multiplying Equation 13 on 
the left by % R?}? and combining the result with Equation 14, 


we obtain: 
w° -W-SRIV (240) 


To convert this into an iterative process, the optimum 
solution is considered to be the value of the weighting vector 


ROtmE Ce Mende ise hole denn 


This equation is referred to as the Newton Method algorithm 
for determining the weighting vector. 

The algorithm takes only one iteration to arrive 
theoretically at the optimum solution for the values of the 
weighting vector. This algorithm would be the ideal method to 
use if exact solutions to the gradient term and the inverse 
correlation of the input vector were known exactly. This is 
not generally the case and an estimate for the gradient must 
be used. To compensate for the lack of knowledge of the 
gradient at iteration k, a constant uw, will be used instead of 
the factor %in Equation 22. The requirements on the range of 


values that p can take on are: 


Ze 


O° ean (23) 


The effect of p on the solution at each iteration is that it 
affects the rate of convergence of the algorithm. A 
discussion of the attributes and further limitations on the 
value of uw will be covered later. A general form of the 


Newton's method of gradient search is therefore: 


W,.. = W,-pR“V, (24) 


This expression is ideal under the following three conditions: 


ei = % 
2. Exact knowledge of the gradient vector, V,. 


3. Exact knowledge of the (unchanging) matrix R°. 


These conditions, unfortunately, are not normally 
attainable. The value of uw is selected between 0 and % so 
that the algorithm 1s overdamped and can accommodate 
fluctuations in the matrix R. The compensation for this 
selection of pp 1s that the algorithm takes longer to obtain 
1€S Optimum solution. A second modification to the value of 
p is to base its value on the average eigenvalue, A.,, so that 
is replaced by pd... 

The second condition pertaining to the exact solution 
to the gradient of the error surface 1s, in general, never 
attainable. The gradient must therefore be estimated in some 


manner. The estimate of the gradient that will be used is 


ee 


based on the square of the present value of the error signal 


and not the true ensemble average as it has been up to now. 


Z 
Oe; = ce (25) 


‘: "ge ae 





From the definition of the error signal: 


aT 
Oe, = Old, BW, = ~X, (26) 
OW 5 OW x 


which leads to the following estimate of the gradient: 


V, = -2e,X, (27) 


With these two relaxations from the three ideal conditzvons, 
the algorithm is now in a form known as the Newton/LMS 


algorithm[Ref 7]: 


Who, = W,+2 pa, R°X, (28) 
This 1S an optimum iterative algorithm that is only complete 
if full knowledge of the input correlation function is known. 
This is generally not the case and thus a final modification 
must be done to produce a working estimate to this algorithm. 
The final modification is an iterative estimate of the inverse 
of the input correlation matrix, R”, and is calleqw@ie. 


Sequential Regression (SER) algorithm. 
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IV. THE SEQUENTIAL REGRESSION (SER) ALGORITHM 


A. DEVELOPMENT 

The Sequential Regression (SER) algorithm uses’ the 
LMS/Newton algorithm as its basis in conjunction with a method 
for iteratively estimating the value of the inverse input 
auto-correlation matrix to produce the weighting vector for 
the filter coefficients. The development of the algorithm 
comes from Adaptive Signal Processing by Widrow and 
Stearns[Ref 7]. The core of the algorithm is the use of an 
estimate for the inverse of the input auto-correlation matrix 
which reduces computational load with minimal estimation 
error. 


The input auto-correlation matrix is 


R=5 (ke x (| (29) 
where the subscript k goes over the entire length of the 
random process X,. R will be estimated iteratively using only 
a finite or truncated sample of the random process X,. The 


estimate of the input auto-correlation matrix is given by: 





at T 
= x. (30) 
R ra . XX; 


This estimate is unbiased when the input variable, X,, is a 
stationary random process. When X, is not stationary, it can 


be shown that the estimate is not very accurate. The estimate 


a 


of R, at each value of k is equally influenced by the 
previous, k-1, estimates of the input auto-correlation matrix. 
In order that the most recent estimates of R, have more 
influence over the new estimate, a "fading" memory term, @, is 
introduced which reduces the significance of the past 
estimates. The fading memory expression of the input auto- 


correlation matrix is denoted by Q, where: 
k 


Q, = oe Gay’, . (38) 


The value of @ is chosen, as a rule of thumb, such that the 
half life of the exponential function is equal to the length 
of stationarity of X,{Ref 7]. This rule of thumb leads to the 


following statements about the value of a. 


a@ = 2(-1/lengthof stationarity of x) (32) 


ee al (33) 


The value of the fading memory factor over k iterations is: 
~a k+l 

yar Syne es (34) 

=0 


and thus the estimate of the input auto-correlation matrix 


becomes: 


k 

= io _ ao k-1 T 

Ro. =, 22222 sop fa (ake ya Ter: (35) 
[oak 1a"? ee = 
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This estimate 1S exact for the condition where X, 1S constant. 
If the input variable, X,, is stationary, @ 1S one, and if the 
limit of Equation 35 is taken as @ approaches one, the result 
1S in agreement with Equation 30. 

The derivation of the SER algorithm begins with the 
estimate of the input auto-correlation matrix. Using this 
estimate in the equation for the optimum weighting factors 


mesules in 
RW, = 68,. (36) 


iijgemaetinition for the estimate of the cross-correlation 


a 


weeeer, P,, 1s, Similarly to Equation 35, 


k 
~ sleet b 4 Pel 
P, = —~—— oe. ox (a7) 
k ae 11 


f=ameenaquations 35 and 37 in 36 gives: 
k 
OM, = Pattd,x,. (38) 
T=0 


The SER algorithm attempts to estimate the next value of 
the weighting vector, W,.,, based on the present values of R, 


and P,. To do this, Equation 38 becomes: 


Ze 


k-1 
QW... = » ga X aed, Xe, 
1=0 


39 
4Q,.,W,+ aX, ome 


(Q,-X, Xi) W, +X ,. 
The last line of Equation 39 is the result of the following 
relationship: 
Q, = “7Q,.,) 32,X;. (40) 

The value of d, can be replaced by recalling Equation 5 and 6 
which state that 

d, = €,+y, = e,+X,W, (41) 
thus Equation 39 becomes 

QW = (Q.-X KW, + (€,+X,W,)k, 


QW, + &%,. 


Multiplying on the left by Q,*, gives an iterative method for 


(42) 


calculating the weighting factor vector 


Wir = Wy + Qee,X,. (43) 
The expression for W,,, in Equation 43 is similar to the 
LMS/Newton method derived earlier. With this in mind, an 
approximation for the LMS/Newton algorithm can be made using 
the above derivation. The SER algorithm used as an 


approximation of the LMS/Newton method is given as: 


hO 
—O 


Die As Se 
Wi = W, + Wi ave | 


1-« Q.&,X,. se) 
The algorithm requires a method of iteratively solving for the 
value of Q.° in order to be complete. 

The method to iteratively solve for Q,' begins by pre- 


feed plving Equation 40 by Q,', post-multiplying by Q,..°, and 


@eetmoy X,, to obtain: 


Qe QQ 1B = Ve GQ QE, + QeKKQ-¥, (45) 


which reduces to, 


Qr1K, = OQ X,+Q, K,XjQ-1%, 


(46) 


Q.X, (a +X1Q;-1%, ) 


The quantity in parentheses is a scalar and can be divided out 


1 


while the quantity X,’Q,.° is also multiplied on the right to 


give 


Qe, 2 kr 


~— = OFX ,X1Qk-1 (47) 
@+X,Q.-1%, 


The first lines of Equation 46 and 47 can be combined and re- 


eenenged to form, 


(Qx-1%,) (Q.°1X,) 4 


= (48) 
% +Xi(Qx-1X,) 


-1 21 
Qk = -=}Qe-17 
Om 


An iterative method of calculating Q,* is now available and is 
only based on its previous value, Q,.,"°, and the present value 


of the input, XK,. The value in parentheses is used three 


zd 


times, it will be calculated separately and represented by 
S.; 

The SER algorithm 1S now completely derived as an 
iterative method of approximating the LMS/Newton method of 
adaptive filtering. The following summary shows the process 
by which each new value of the weighting factor vector is 


calculated. 


a2 (-1/lengthof stationarity of X) 


e. = d,-W,X; 


-1 1 ee a 
Qx = q(Qe1- 887] 


Ze 
ea cn gie,X, 
1 
O< p< i UA eee 





max 


The values that must be carried forward from one iteration to 
the next are Q,.,°, and W,. These represent the previous value 
for the estimate of the inverse input auto-correlation matrix 
and the present estimate for the weighting factor vector. 
The output of the algorithm is the error, €,, between the 
present value of the desired signal and estimate of the 


desired signal: 
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OUTPUT = €, = d,=W,X; (50) 
This value is saved after each iteration in a vector 
representing the uncorrelated portion of the desired signal 
and noise vector of the primary channel with the noise only 


vectors of the reference channels at the time index k, 


X,.. [Ref 7] 


B. INITIALIZATION 

The SER algorithm assumes that the value of the inverse 
input auto-correlation matrix and the weight factors vector is 
known from the previous step. At some point in the past, an 
estimate for these terms must have been made so that a 
starting point can be established. The required accuracy of 
the estimates depends on the length of time that the algorithm 
will be operating and the earliest time that accurate data is 
required. 

The driving force behind the accuracy of the initial 
estimate is the speed with which the algorithm can converge to 
an optimal solution. The speed of convergence is dependent on 
two elements of the algorithm. The first element is the 
initial value of the inverse input auto-correlation matrix, 
Q°. The significance of this initial element is that the 
further the initial estimate is from the true value, the 
further the algorithm has to adapt in order to approximate the 


optimum solution. The second element is the step size, 2upA..., 


Bd 


which specifies how far the algorithm can progress towards the 
optimum value in one iteration. The step size will be 
discussed in the next section but for now it shall be assumed 
to be a constant in this discussion of the initial conditions. 

The use of a large constant multiplied by the identity 
matrix as the initial value of the 0, has been sufficienely 


argued by Lee [Ref. 8]. 


Q5° = QI (51) 

As discussed by Lee, the use of a large constant for q, 1s the 
proper selection for the case where the random process is 
Stationary. Widrow and Stearns argue that this same choice 
for the initial value of q, in anon-stationary condition, is 
also appropriate[Ref 7]. The point is made that if a priori 
knowledge of the random process allows for a better estimate 
of the value of Q,', then that information should be used. 
The method of initialization for this study is to use Eova@mien 
51 with q, equivalent to 100. 

The initial values used in the weighting factor vector are 
assumed to be zero. This assumption is made since knowledge 
of the process before time zero is unknown. The use of a 32 
tap delay filter also allows for an estimate of the initial 
full weight factor vector but this is done so at the expense 
of the first 32 output data points. This is done by expanding 
the number of elements in the weight factor vector as each of 


the first 32 points are read into the algorithm. The 


algorithm 1s applied only to the data that is available; 
therefore, the first time that each element of W, is non-zero 
is at time index 31. This corresponds to the time when the 
first 32 data points are input into the algorithm. The result 
of this initialization is that the first real data point to be 


considered an output of the SER method used here is at k = 31. 
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V. STABILITY AND CONVERGENCE 


A. THE SINGLE WEIGHT CASE 

The SER algorithm contains values that can be set to 
affect the rate at which the algorithm achieves the optimum 
solution. A simple case where only one weight factor is 
concerned will be used to demonstrate the relationships 
involved in the selection of the convergence rate factors. 
The single weight iterative process using a gradient search 


method like the SER algorithm can be represented by: 
Wiad = Wet Hh (-Vy) (52) 


The expression for the single weight gradient at time 


index k is given by: 


a 
Vi. - 2 
a w=W, 


. GdlA(w-w?)?] 
dw 


(S37 





w= WwW) 
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The value of the input auto-correlation matrix, Yp,. for the 
Single weight case, is replaced by its equivalent eigenvalue, 
A=E{x?]. Combining Equations 52 and 53, an iterative eqmeammem 


for the weight factors is found: 


Wii = W,- 2pA (w,-Ww?) (54) 
for the single weight case. The terms in Equation 54 can be 


rearranged to combine like terms. 

W,=w°t (1-2pA) *(w.-w?) (55) 
This expression contains the geometric ratio of successive 
terms which is also used to define the stability requirements. 


Stability 1S guaranteed if the absolute value of the 


geometric ratio is less than one, i1.e., 
peo <1 (56) 


The limits placed on the value of the uw can be expressed in 
terms of the eigenvalue of the input auto-correlation matrix 


A by rearranging Equation 56 


E 


FH 0 (oy) 


The optimum solution can be seen as the solution of Equation 
55 if the limit of w, as k goes to infinity. This 1s because 
the geometric ratio approaches zero as k approaches infinity. 
The rate at which the algorithm approaches the optimal value 
1s dependent on the value of the geometric ratio. 

Figure 4 shows the effect of varying values of the 
geometric ratio on the rate of convergence. If the value of 
the geometric ratio can be maintained between zero and one, 
then the algorithm will always approach the optimal solution 


from one side and can be considered as overdamped. EOi? ia 


a0 


geometric ratio of zero, the algorithm reaches the optimum 
value in one iteration from one side and can be considered as 
critically damped. If the magnitude of the geometric ratio 
falls between zero and negative one then the algorithm is 


underdamped. 
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Figure 4 Effect of Geometric Ratio (g) on the Rate of 
Convergence: overdamped O<g<1; critically damped g=0; 
underdamped -l<q<0 
B. MULTI-WEIGHT CONDITION 

The single-weight condition discussed in the previous 
section can be expanded to develop an understanding of the 


multiple weight factor case. The same conditions that applied 


to the single weight case are true in the multiple weight 


oro 


@eee The difference is that each weighting factor has an 
associated eigenvalue and thus the condition of stability is 
now dependent on a set of simultaneous equations that must 
meet the stability condition. 

(1-2pA,) 

(T=2N Aa) 


(1-2A,) 


W=Wor (W,-W°) (58) 


(1-2p4,_,)| 


The subscript on A indicates the index value corresponding to 
a specific weight value. There are a total of L weighting 
factors and thus there are L individual eigenvalues, A. 

The value of p in each of the simultaneous equations is 
held constant. Since u is the same for all equations, the 
eigenvalue that has the largest value will impose the most 
restrictive conditions on the selection of a value for yp. The 
new range of values that uw can take on is: 


1 


max 


>p>0 (29) 
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The condition stated in Equation 59 gives the multi-weight 
Stability requirements that must be met to ensure that the 
algorithm will converge to the optimal solution. 

The rate at which the algorithm will converge was shown in 
Figure 4 to be dependent on the value of the geometric ratio. 
Now that the value of wu is a function of the maximum 


elgenvalue of the input auto-correlation matrix, R, it can be 


of 


shown how this affects the rate of convergence. By holding 
the value of wp constant for each weighting factor, an 
undesirably slow convergence rate can occur for widely varying 
values of A. As shown previously, the algorithm converges in 
one iteration when p is equal to 1/2A for the single weight 
case. Using that same criteria, the value of the multi-weight 
constant p is 1/2A,... The use of this criterion ensures that 
the algorithm falls into the overdamped category since yp is 
always less than or equal to 1/2A, for 1=[{0, 1, 2, J) 3aie 
The weight factor, w,, with the smallest value of A will have 
the slowest convergence rate since its geometric ratio will be 
farthest from zero. This is true because the value of 2pA 
will be at itsS minimum and will cause the geometric ratio to 
approach one. When this condition holds, the algorithm may 
take excessively long periods of time to converge to an 


optimal solution. 


C. NON-STATIONARITY EFFECTS 

The discussion up to this point in dealing with stajbuie 
and convergence assumed that the process under study was 
stationary. Under stationary conditions, the values of the 
eigenvalues of R are assumed to be constant. A second key 
assumption to the development is that the values of the 
elgenvalues are known. The lack of knowledge of the input 
auto-correlation matrix, R, requires that some estimation be 


made for the value of p. The value of pp must be adjusted to 
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ensure convergence of the algorithm in an environment where 
the eigenvalues of the auto-correlation matrix are varying at 
some unknown rate from iteration to iteration. 

A method for compensating for the variations in the range 
of values that p can assume 1S to use an averaged eigenvalue, 
Awe, Lor the selection of the constant w for each iteration. 
The value of p 1S now a constant for a single iteration but is 
changed from iteration to iteration based on the average 
eigenvalue. Stability of the algorithm is guaranteed if the 
value of the product, wiA,,., 1s between zero and one. It is 
important to note that the expression pA, 1s now used only to 
replace the constant p used previously. The idea is to scale 
the average eigenvalue for that iteration by a constant, un, 
and use it as the factor for controlling convergence of the 
eeegocit hm. 

The geometric ratio which controls the convergence of the 
algorithm has now been altered to include a term to compensate 
for the non-stationarity of the input random process, X. The 


new expression for the geometric ratio is 


ieee Ret (60) 


ave 
imemmmverse input auto-correlation matrix, R’, is included in 
the development of the SER algorithm. The convergence 
controlling term is thus nwpd,,., and must be less than the 
maximum eigenvalue of R to ensure stability and convergence. 


A second condition that can be imposed on the product pA... 1s 


oo 


that its value must always be between zero and one. Using the 
fact that the spread between maximum and minimum eigenvalues 
can be rather large, the product of wpdA,,. must be much less 
than one. 

The ELF background noise processes have shown a tendency 
to have a wide range between eigenvalues and to slowly change 
over time. Because of these two conditions, the value of the 
Aw. Must be adjusted over time to account for these 
Variations. The chosen method is to estimate the average 
elgenvalue by taking the average of the square of the input 
VECl Olas. This value is multiplied by the 1/100 of the 


inverse of the maximum of the squared input vector. 


S x? 
1 1=0 


: (61) 
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The use of this estimate for convergence control in the 
algorithm ensures that the requirements for stability are met 
for each iteration. The 1/100 term ensures that the maximum 
value of p 1s not exceeded. It was derived empirically with 
the use of an actual set of ELF background measurements and 
may require adjustment ina different environment. 

The choice of using the present input to the algorithm to 
get an estimate of the average eigenvalue is used for 
computational ease with minimal added error. The use of the 
inverse input auto-correlation matrix should provide a more 


accurate measure of the eigenvalues; That, however, requires 
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the inversion of the matrix. The prevention of matrix 
inversion is the exact reason the SER algorithm is developed, 
so to obtain the eigenvalues in this way defeats the reason 
for the algorithm's design, and is therefore counter- 
Precuctive. The method has proven adequate and reasonably 
accurate in empirical work up to this point, but does leave 


room for improvement and further study. 
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VI. IMPLEMENTATION 


A. A PHYSICAL DESCRIPTION 

The Applied Physics Lab has two identical sensor platforms 
submerged in shallow water. Each platform has electrode pairs 
for detecting the geoelectric ELF background noise. The 
system has been in operation for over two years collecting 
samples of the background noise environment. 

Ideally, the two platforms are parallel when the platforms 
are placed on the bottom of the bay. Figure 5 shows the ideal 
alignment of the platforms. The two platforms are labeled as 
the east, E, and west, W, platforms. Each platform has a set 
of orthogonal electrode pairs, labeled X and Y. These 
electrode pairs are orientated to receive the two orthogonal 
source vectors, E, and By | in ene wnorizoneatee sane The 
strongest geoelectric field strengths are in the horizontal 
plane due to the shift in polarization as a wave crosses the 
water-alir boundary eeee vertical to horizontal. 

If the two platforms are perfectly aligned and the noise 
recorded on each are identical, the noise in a specific 
direction, X or Y, can be removed by subtraction. To cancel 
the x-direction noise simply subtract XW from XE. Because the 
platforms are not perfectly parallel and the recorded noises 


are not identical though highly correlated, this method 1s not 
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Figure 5 Illustration of Sensor Set-Up 


ideal. The actual platforms can be misaligned up to 10 
degrees from perfectly parallel. This is due to physical 
constraints such as bottom contours and mounting technique. 
The result of such a misalignment is that both the XW and YW 
sensors of the west platform have components that are coherent 


with XE of the east platform. If the electrode pairs in 
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Figure 5 are thought of as vectors, then it is easily seen how 
misalignment can cause coherency between XE, XW, and YW. 

A set of real data is traced through the signal processing 
system designed to accompany this sensor system. A calibrated 
source signal is embedded in the signal XE. The path that the 


received signal follows can be seen in Figure 6. The output 
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Figure 6 The ELF Signal Processing Block Diagram 
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of each sensor 1s sampled at 128 Hz. The data is bandpass 
filtered and resampled down to 32 Hz for processing. The data 
is then passed through the SER algorithm where the coherent 
portions of the west platform are removed from XE. The output 
of the SER portion of the processor is demodulated and passed 
through a likelihood ratio detector where a decision can be 


made as to when the embedded signal is present. 


B. INPUTS TO THE SER ALGORITHM 

The SER algorithm is used to reduce the ELF background 
noise detected in the primary submerged electric field sensor 
by cancellation with the output of the reference parallel 
electric field sensor. These parallel reference sensors are 
also submerged and are, ideally, assumed to be located within 
the same ELF background noise as the primary channel. The 
primary channel y represents sensor XE and the reference 
channels, xl and x2, represent XW and YW, respectively. The 
samples are the actual output of the sensor system used by 
APL. The samples are taken at the same time without any known 
Signals present. The sampling rate of the data input for the 
SER algorithm is 32 Hz and the total length of data is 113760 
points. These data samples correspond to almost an hour of 
background noise measurements. 

The filter is a 32 point transversal filter that is 
applied to each channel simultaneously. The length of the 


transversal filter means that the previous 31 samples of the 
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sensor are used with the present sample value for each 
1teration of the algorithm. The goal of the algorithm is to 
take the output of the two reference channels, xl and x2, and 
to generate an estimate of the noise in the primary channel y. 
Figure 7 shows what the primary channel's data looks like in 
real time. The length of stationarity, used to determine the 


forgetting factor, @, is assumed to be about 300 data points 
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Figure 7 The Primary Noise (Y) and Calibrated Source Signal 


long. This value is used based on empirical analysis of the 
real data. 
The effectiveness of the algorithm is analyzed by running 


a set of scenarios with a known signal inserting in the data 


Popes. The scenarios consist of manually inserting a 32 Hz 
Siemal On a 5S Hz carrier that is 8192 data points long. 
Figure 7 shows the components of the primary channel signal in 
a real time representation. The primary channel is processed 
by the SER algorithm while varying the strength and location 


of signals embedded in the background noise y. 


C. SIGNAL STRENGTH RESULTS 

The algorithm is first run with the location of the signal 
within the noise held constant. The strength of the signal is 
varied by dividing each data point in the signal vector by a 
Censcant. A series of four runs are conducted, each 
succeeding run has half the signal strength of the preceding 
run. The initial run uses the signal shown in Figure 7 with 
the amplitude divided by 500. The signal is indistinguishable 
prior to processing as seen in Figure 8. 

The signals for each run are centered at time 19.2, 32 and 
44.8 minutes into the noise sample sy. There 1s no 
Significance to the positions chosen except that the signal 
data streams were not allowed to overlap. The objective of 
this first set of runs is to demonstrate how the SER algorithm 
enhances the ability to distinguish small signals buried in 
the noise. The data is first processed in the SER algorithm 
to remove the coherency between the reference channels and the 
primary channel. The output is then demodulated leaving only 


the in-phase and quadrature components of the non-coherent 
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portion of the primary channel. An identical reference set of 
data is demodulated without passing it through the SER 
algorithm so that a visual and quantitative analysis of the 
effects of the SER algorithm can be made. The output of the 
demodulator is shown in Figures 9 and 10. Note that the 
envelope of the signal is visible to nearly the strength of 
1/2000, ov the tha rdvrun) made 

The output of the demodulator is passed as input into a 
matched filter detector based on the demodulated signal. The 
output of this detector 1s shown in Figures 11 and 12 within 


SER and non-SER output for a given signal strength shown 
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Figure 8 Input Primary Channel Noise And Signal 
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together. Numerically, the peak value out of the detector is 
three to five times greater when the SER algorithm is used as 
compared to when it is not used. 

A final numerical analysis is performed by uSing a primary 
channel without a Signal installed to get a measure of the 
noise present in the output. The standard deviation of this 
output of the detector when this noise only vector is used as 
input 1s used aS an approximation to the noise power. 
Dividing the outputs of the detector in Figures 11 and 12 by 
this estimate of the noise power, an estimate of the signal- 
to-noise-ratio (SNR) can be made. The values of the average 
peak SNR for the three locations in each run are plotted 
versus Signal strength in Figure 13. The SNR of the SER 
applied data runs about 6 dB higher than the non-SER applied 


data. 


49 


Amplitude 


—0.2 = 


O LO” eet Z0 uD) ice) 


Time (min) 


Amplitude 


O Oe G26. 30m > SO 
Time (min) 
b 


Amplitude 


O 20° 30 40° 5 


Time (min) 


Amplitude 


O 10.20" 2350) ae 
Time (min) 
d 


Figure 9 Demodulator Output of 
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Figure 10 Demodulator Output of Quadrature Components; 
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Figure 11 In-Phase Detector Outputs: 
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Figure 12 Quadrature Detector Outputs: (a) Signal/500 
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Figure 13 SNR vs Relative Signal Strength 


D. DETECTOR THEORY 

The process of adaptive noise cancellation is performed 
for the purpose of increasing the signal-to-noise ratio. A 
section on the narrowband detection schemes that have been 
considered in conjunction with the noise canceler are 
presented since they are integral to the analysis here. Two 
likelihood ratio detectors are considered for the detection of 


narrowband signals. 


The first detector is designed with the assumption that 
not only is the signal known, but its peak within a given 
interval is constant relative to the start of that interval. 
This last assumption 1S to assure that the demodulated in- 
phase and quadrature signals that are used in the likelihood 
ratio detector match those in the noise canceled and 
demodulated output of the sensor system described previously. 
In the second detector, the latter assumption was dropped and 
the design followed that of an optimal envelope detector. [In 
both cases it was assumed that the residual noise field was 
gausSian. This assumption, while not strictly true, seems to 
yield quite reasonable results.[{Ref. 9] 

The signal is generated by a calibrated source with known 
qualities. The modelling equation for the narrowband signal 


1S given by: 
Suet) Cos | 2 eet +104 t)] (62) 


For the first detector, the phase term, 0, is assumed constant 
within any non-overlapping interval 0 to T. The interval of 
0 to T 1s equivalent to the length of the signal. The purpose 
of this assumption is so that a Matched Filter can be created 
for the in-phase and quadrature components by using the 
demodulated in-phase and quadrature component of the signal. 
Figure 14 shows the demodulated signals used for the detector. 

The first detector is a multi-channel Matched Filter shown 


mm gure 15. The Matched bilters, fh, and h,, are determined 
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Figure 14 Demodulated Signal for Match Filtering 


based on the requirement that the output signal-to-noise ratio 


1S maximized at time T. The corresponding equations are 


7 [Ri (t-t) Ry, (t-t) |[A, (+) 
Jlrs (t-2) ee (C20) qe A) 


where the quantities R denote the various auto and cross 


S. (iw) 
ea (63) 








correlations. This detector is the one used for all 
processing performed in this thesis. The input to the 
detector is the noise cancelled output of the SER portion of 
the processor after demodulation. [Ref. 10] 

The second detector iS an optimal nonlinear envelope 


detector. The description of its principles and design are 
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Figure 15 Linear Multi-Channel Match Filter 


located in the Appendix. The detector was not used by the 
author in this work relating to the analysis of the SER 
algorithm, but it 1s presented as complementary work relating 


to the overall project. 


E. SIGNAL LOCATION 

The signal strength analysis showed that the signal 
divided by 2000 is easily discernable from the background at 
the output of the detector. The effect of varying the 
location of the signal will allow analysis to get a better 
feel for the expected SNR based on a larger sample size for 
averaging than previously used. A second benefit of this 
analysis is to see the effect of the initialization process on 
the SNR. The signal will be shifted by a quarter of its 
length for each run, or 2048 points. The SNR is plotted 
versus the location of the center of the signal in Figure 16. 

The data is run with and without the SER algorithm to get 


a good feel for the enhancement achieved with the SER 
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algorithm. Three separate locations of the signal will be 
used for each run to reduce calculation time. The use of 
multiple signals causes no effect on the SNR calculation as 


long as the signals are not allowed to overlap. 


VII. RESULTS AND CONCLUSION 


A. ANALYSIS OF RESULTS 

The use of the SER algorithm enhances the ability to 
detect narrowband low power signals buried in the background 
noise. Figure 17 shows the spectrum of the noise sample from 
sensor XE before and after the use of the SER algorithm. The 
frequency range is from one to ten hertz, which is in the 
range of concern for the JHU project. The spike seen at four 
hertz (mote the x-axis 1s a logarithmic scale) 1s caused by 
the processing equipment. This 1S a man-made noise scurce and 
will be filtered out in the final project design. The figure 
illustrates the 10-15 dB of noise cancellation possible by 
using this method of adaptive noise cancellation. 

The location of the signal in comparison to the start of 
the adaptive process has been shown to have very little effect 
on the signal enhancement. The system implementation will 
most likely cause the algorithm to only require initialization 
once a week and operate constantly otherwise. The 
initialization process takes about 100 iterations or 1.67 
minutes before settling out but is heavily data dependent. 
Even while settling, the performance of the algorithm is 
adequate to receive most signals. For the purpose of this 


M@eenect, initialization should be of no concern. 
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Figure 17 Illustration of SER Algorithm Cancellation 
Capability, Power Spectral Density of 1-10 Hz 


B. CONCLUSIONS AND RECOMMENDATIONS 

It was hoped at the beginning of the research that 
information could be gained on the characteristics of the 
noise environment. The stationarity of the background noise 
environment was found to be at least nine seconds. This is 
believed to be only a characteristic of the specific data set 
that 1s used. The length of stationarity affects the value of 
the forgetting factor, a, and 1s implied by the use of the SER 
algorithm. A better feel for a good average value to use for 


the length of stationarity will result from the use of more 


60 


data sets. This analysis needs to be conducted and could be 
a source of further work in this area. 

A second area of work that could develop from this 
research is a study of the cancellation possible from a 
vertical above-surface sensor and a submerged horizontal 
sensor. The purpose of this would be to differentiate between 
surface and submerged activity. This may prove to be 
impossible but it would be interesting to see the amount of 
correlation possible between the two locations. This type of 
arrangement could also lead to a wider distance between 
sensors which could also increase the knowledge of the 
background environment. 

The purpose of the research presented in this thesis is to 
develop a method to adaptively cancel the ELF noise present in 
a submerged sensor array. The SER algorithm appears to be an 
adequate method to accomplish this goal of noise cancellation. 
The cancellation ratio in the frequency range of interest is, 
on average, 12.5 dB. This cancellation ratio along with the 
use of optimal detectors will assist the system in detecting 


narrowband signals buried deep in the background noise. 
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APPENDIX 

The second detector used in conjunction with this project 
of measuring the geoelectric background noise is an optimal 
nonlinear envelope detector. The symbolism used is slightly 
different but efforts are made to define all relevant terms. 
The in-phase and quadrature data (after noise cancellation) 
are denoted by rj and rf, respectively. Similarly the in-phase 
and quadrature noises are represented by ni and nf. Assuming 


a joint Gaussian distribution we have the following two 


hypothesis: 
Hy 2 ( 2 QO) eer 
Hi os € ro, )) =n? Sen ee ee ec 


where H, denotes the null hypothesis (noise only) and H, 
denotes the alternative (signal and noise present). The 


vector notation means 


Ze) eer er ere 
where N denotes the length of the interval for which the 
decision is being made. It is then clear that the probability 


density function for the received data is as follows: 
a ee eer qo a 
p (r|H,) exp| = (r-s)"R> (r-s) 


where R,, denotes the covariance matrix of the noise vector. 


In this particular case, it 1s assumed that the in-phase and 
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the quadrature components have been decorrelated and whitened 
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the above over the angular variables will give 








N 
ty N N 
Z i ee SI, 
ip CEH.) =we 2 exps= (pees .) a 
Ms oma 0? 2 a I] ‘\ 0? 
where I, 1s the modified bessel function of order zero. The 


vector ron the left hand side now only refers to the various 
envelopes, r= [fnr,, ..., ©]. The same is true for the signal 
envelope. The likelihood ratio detector is found to be: 

Rp Cre, ) 


= =eXD 
P(zr|H,) 











aac : Si2, 
20-4 = G* 


The log-likelihood ratio is then: 





N N 
SG 
loga=-— de &x* log rf x : 
20° =7 =1 o* 
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The first term is a constant, dependent only on the known 
Signal, and can be dropped. Finally the test statistic or 


detector output is[Ref 10]: 





= Sr 
LOgA= s LOgnWi ae 
dy { o? 
The implementation and analysis of this detector scheme is 


left as a possible follow-on thesis topic. 
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